
3 


76 



00163 



/\jf\-sn 77Tl'<5>/35t) 


NASA-TM-81350 19800015607 


NASA Technical Memorandum 81350 

1 

A 


BIVARIATE NORMAL CONDITIONAL AND RECTANGULAR PROBABILITIES: 
A COMPUTER PROGRAM WITH APPLICATIONS 


Ram Swaroop, James D. Brownlow , George R. Ashworth and William R. Winter 


May 1980 


c 


J 


! f r; ,P j 


in 


MAY ii ii 


U--r- 


Liti: j ; 





NASA Technical Memorandum 81350 


BIVARIATE NORMAL CONDITIONAL AND RECTANGULAR PROBABILITIES: 
A COMPUTER PROGRAM WITH APPLICATIONS 


Ram Swaroop , James D . Brownlow and George R . Ashworth 
System Development Corporation 
Edwards , California 

and 

William R. Winter 

NASA Dryden Flight Research Center 
Edwards , California 


tvi/\s/v 

National Aeronautics and 
Space Administration 


1980 


I0K(j-3-‘^09 ^ 54^ 



INTRODUCTION 


In applications involvinq univariate data where estimates and confidence 
intervals are required, the normal distribution is commonly employed. This 
distribution is mainly utilized because the probabilities under a normal 
curve are readily available. In contrast, use of multivariate probabilities 
in p-variate normal data are less frequent, primarily because probabilities 
for the multivariate normal case are generally not available. Except for 
very special cases, the probabilities for sections of p-dimensional space 
require extensive computations, since the canonical multivariate normal den- 
sity changes with every change in correlation coefficient parameters. Even 
the probability computation in the bivariate normal case (p = 2) with only 
one value for the correlation coefficient over arbitrary sections of the 
(x, y) plane is not easy. Probability computations, therefore, in p > 2 
dimensions are correspondingly much more difficult. (Ref. 1) 

In many applications, problems are posed which not only require the^ 
probabilities over a section of p-dimensional space, but also the conditional 
probabilities of r (r < p) variables when the remaining (p - r) variables are 
either fixed, or are within designated intervals. For example, in aircraft 
target tracking studies, it is of interest to know the probability of X 
deviations from the target when Y deviations are considered within desig- 
nated bounds. In aircraft performance studies it is important to know the 
distribution of the pilot's cardiac R-R intervals either under an assigned 
difficult aircraft maneuver or under the dynamic flight conditions. 

The results on conditional and marginal distributions of r variables 
when the (p - r) remaining variables assume fixed values are well estab- 
lished. (Ref. 1) Similar results, when the remaining (p - r) variables 
assume values within specified ranges involve complexities and are discussed 
in this report. 

In this study, results on bivariate normal distributions (p = 2) are 
reviewed. Various derivations and properties of bivariate normal conditional 
probabilities are derived. A computer program for conditional probabilities 
for all assigned values is included. From conditional and marginal proba- 
bilities, the rectangle probabilities are then obtained. Examples are pre- 
sented to illustrate the use of the program. The program listing is appended 
to this report. 


SYMBOLS 

Ay lateral acceleration 

A^ vertical acceleration 

c a constant with fixed numerical value 

exp(x) exponential function at x 



F(s) 

f(u, v) 
f(x). f(y) 
f(x, y) 
f(x]a<Y<b) 

f(x|Y=y) 
f(x|Y<t, p<0) 

f(x|Y>-t, p>0) 

Gp(s. t) 

9t(x) 
g^(x|p>0) 
g^(x|p<o) 
p. r 

Pr[a<Y<b] 
Pr[c<X<d, a<Y<b] 

Pr[X<h, Y<k] 

U, V, X, Y 
u, V, X, y, t 

\ 

a 


conditional distribution of X at X = s Given Y is 
in interval (a, b) 

general bivariate normal density 

standard normal densities 

standard bivariate normal density at X = x, Y = y 

conditional density of X at X = x given Y is in 

interval (a, b) 

conditional density of X at X = x given Y = y 

conditional density of X at X = x given Y is less 

than t and correlation is negative 

conditional density of X at X = x given Y is greater 

than -t and correlation coefficient p is positive 

double integral with two arguments s and t vnth a 
fixed value of correlation coefficient p 

conditional density of X at X = x given Y is in 

interval (-t, t) 

conditional density of X at X = x when correlation 

coefficient p is positive and Y is in interval (-t, t) 

conditional density of X at X = x when correlation 

coefficient p is negative and Y is in interval (-t, t) 

dimension of multivariate data or distribution 

probability that variable Y is in interval (a, b) 

joint probability that variable X is in interval (c, d) 
and variable Y is in interval (a, b) 

probability that X is less than h and Y is less 
than k 

random variables 

specific values of random variables 
forward velocity 

fixed positive constant less than 1 
mean of forward velocity 
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a 


c 


'u’ 


"y 

*(t) 


mean of subscripted random variable 
mean of lateral acceleration A 

J' 

mean of vertical acceleration 

correlation coefficient between two random variables 

standard deviation of forward velocity 

standard deviation of subscripted random variable 

standard deviation of lateral acceleration Ay 

standard deviation of vertical acceleration A^ 

standard normal distribution at t 


BIVARIATE NORMAL DISTRIBUTION 

A bivariate normal distribution of a random vector (U, V) is charac- 
terized by parameters: and p. The density function 



is defined over the entire (u, v) plane. When the variables U and V are 
standardized, by defining the new variables 

f(x, y) = ^2 tt/i - exp|^- (x^ - 2pxy + y^) / 2(1 - p^)J 

the density function of (X, Y) reduces to the canonical bivariate normal 
density 

X = (U - , Y = (V - 

defined over the entire (x, y) plane. The parameter p is called a cor- 
relation coefficient and takes values in the interval (-1, 1). Without any 
loss of generality, this canonical density f(s, y) is considered in this 
study. 
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The density function f(x, y) exhibits certain properties. It is 
symmetric in opposite quadrants since 

f(x, y) = f(-x, -y) 

and 

f(x, -y) = f(-x, y) 

Further, f(x, y) is constant over all the ellipses 

x^ - 2pxy + y^ = c(l - p^) 

for every value of x. (Fig. 1) The intercepts made by these ellipses on 
the X and y axes are equal. If p is positive, the major axis of the 

ellipse is along the 45“ line with a length of 2/c(l + p); and the minor axis 

is along the 135“ line with a length of 2/c(l - p). If p is negative, the 

major axis is along the 135“ line with a length of 2/c(l - p); the minor axis 

along the 45“ line has a length of 2/c(l + p). (Ref. 2) The ellipse 

x2 - 2pxy + y2 = (1 - p2) log 1/(1 - a)2 

for all 0 < a < 1, contains the a proportion of the (X, Y) distribution. 
(Ref. 3) 

The marginal distributions of X and Y are standard normal with the 
covariance between x and y equal to p. When p = 0, then 


f(x, y) = (/2 tt)“^ exp( -x^/2) (v^)"^ exp( -y^/2) 

= f(x) • f(y) 

which is a product of standard normal densities, implying that p = 0 if and 
only if X and Y are independent. When p f 0, bivariate normal probabil- 
ities Pr(X < h, Y < k) for a few selected values of h and k are avail- 
able from tables and graphs. (Ref. 4, 5) For general values of h and k 
approximation and interpolation methods are used. 

DERIVATION OF CONDITIONAL DENSITIES 

Conditional Dengitv nf X Gi ven Y = y . It was stated earlier that if a 

random vector (X, Y) has a bivariate normal distribution, then the marginal 
distribution of either X or Y is normal with mean 0 and variance 1. 

The conditional distribution of X for a fixed value of Y = y, however, is 

normal with mean py and variance (1 - p^). The conditional density f(x|Y = y) 
is derived below. 
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Figure 1. Marginal, conditional densities and ellipses of constant 
densities from bivariate normal density . 


f(xlY = y) 


which is the 
(1 - p^) and 

Conditional 
given a < Y 

f(x|a < Y < b 


= f(x. y)/f(y) 

_ ( 211/1 - p^) exp[-(x^ - 2pxy + y^) / 2(1 - p^)1 

exp( -y^/2) 

= [/2it( 1 - p2)J exp |-[x^ - 2pxy + y^ - (1 - p^)y^J j 

= [/2it( 1 - p2) ] exp[-(x^ - 2pxy + p2y2)/2(l - p2)] 

= \/zAl - exp[-(x - py)V2(l - P^)] 

I 

density of a normal distribution with mean py and variance 
is shown in Figure 1. 


ipnQitv of y Rivon a < Y < b . The conditional density of X 
< b is not normal and is derived as follows. 


, _ (-^r r exp |'-(x^ - 2pxy + y2)/2(l - p2)] dy 
(v^) exp( -y^/2) dy 

= [$(b) - $(a) ] ^ ^2ir/l - ^ • 

/ b 

exp - 2pxy + p^x^ + x^(l - p2)] ^ 2(1 - p2) I dy 

= exp{ x^/ 2 ) [$(b) - $(a)] 

/[^ exp[-(y - px)2/2(l - p2)J dy 

= f(x)[$(b) - ${a)]'^ |$[(b - px )/ A - p^] - $[(a - px)//l - p2]| 
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where 


f(x) = (v/2¥)"^ exp( -xV2) 


is a standard normal density and 



is the standard normal distribution function. 

This conditional density is neither normal, nor symmetric. However, in 
special cases discussed below, symmetry is identifiable. 

Civmmptrv in r.nnditioninq -t < Y < t . With -t < Y < t, the conditional 
density of X at specific values of x and -x are 


9t(x) 


f(xl-t < Y < t) 

f(x)[$(t) - (-t)] U(t - px)//l - 


'[(-t - px)//l - P^] I 


9^(-x) = f(-x|-t < Y < t) 

= f(-x) [4>{t) - ${-t)] + px)//l - p^] - $[(-t + px)//l - p^] I 


The symmetry of a standard normal density shows that f(“X) ~ f(x). With the 
asymmetry of distribution function $(t) = 1 - 4>(-t), it is seen that 


$[(t + px)/ A - p^] - <i>[(-t + px)/ A - p^] 

= 1 - 'I>j^(-t - px) / A - - |l - $[(t - px)/ 

= $[(t - px)/ A - p^] - $[(-t - px)/ A - p^] 


- p 




Thus g^(-x) = g^(x), showing that for -t < Y < t the conditional density 

of X is symmetric in x, as shown in figure 2. 

The conditioning, -t < Y < t, with positive and negative values of 
correlation coefficient p also show symmetry of g^(x). It is to be noted 


that 


g^(x|p > 0) = f(x)[<i>(t) - $(-t)]“^ 


[(t - px)/ /l - p^_ - $[(-t - 


px 


)/.TT7 
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Figure 2 . Conditional density of X given -t < Y < t (t = 1 . 000 , 
probability = 0.6826) where (X, Y) is bivariate normal with p = 0.9000, 
and standard normal density . 




g^(x|p < 0) = f(x)[$(t) - $(-t)] ^ 



?] - $[(-t + px)/ A - p^]| 


By the symmetry of f(x), the asymmetry of $(t), and the arguments given 
earlier, it is seen that g^(xlp > 0) = g^(x|p < 0). The graph of such a 

density is shown in figure 2. 

<;vmmgjtrv whf»n -oo < Y < t and -t < Y < +« . In these cases it is to be noted 
that »(-“) = 0, $(«) = 1. Thus the conditional densities of X are 

g^(xlp > 0) = f(xlY < t) 


= f(x)[$(t)] 



px 


)/ A - p^] 


g_^(x|p > 0) = f(x|-t < Y) 


= f(x)[l - $(-t)] 

= f(x)[$(t)] - P^] I 



g_^(-x|p > 0) = f(-xl-t < Y) 

= f(x)[$(t)] |'^’[(t - px)/ A - p^] I 

Thus g^(x) = g_^{-x), showing that a one-sided conditioning on Y yields 

the same density for x as does the conditioning on the other side for the 
opposite X. Further, for negative and positive values of p, it is to be 
noted that 


g^(x|p > 0) = f(x)[<i>(t) ] |<j>[(t - px)/ A - p^]| 

and 

g^(x|p < 0) = f(x)[$(t)] |$[(t + px)/ /l - p^] I 


= g_t(x|p > 0) 
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Therefore, if the conditioning on Y and the sign of the correlation co- 
efficient are reversed, the density remains invariant. An example of these 
densities is shown in figure 3. 

DERIVATION OF CONDITIONAL DISTRIBUTIONS 

Conditional Distribution Function nf X Rivon Y = y . The distribution 
function from the conditional density 

f(x|Y = y) = [/z (1 - ] exp[-(x - py)^/2(l - p^)j 

derived earlier, is easily obtainable via the normal distribution function 

with mean py and variance (1 - p2). it is to be observed from figure 1, 
that mean py is a function of the correlation p and the specific condi- 
tioned value of y, but the variance depends only on p and is invariant 
for all values of y. Thus the width of any a level confidence interval 
remains the same irrespective of the conditioned values of y. 

In applications, the conditioning of variable Y is seldom a fixed 
value. The conditioning is usually in a range a < Y < b, and the formulae 
for this case are different from the results for Y = y. 

Conditional Distribution of X Given a < Y < b . The conditional density 

f(x|a < Y < b) = f(x)[$(b) - <i>(a)] j$[(b) - px)/ A - p^] - <f[(a - px)//l - p^]| 

where 

f(x) = (v'2tT)"^ exp( -x^l2) 
and 


$(t) = / f(x)dx 

*/ —00 


was derived earlier. A general expression for the distribution function 

F(s) =J^f{x\a < Y < b)dx 
-00 

= [$(b) - $(a)]“^f f{x) j $[(b - px)//l - P^] - o[(a - px)//l - p^] | 
^00 


dx 


for all the values of s involves integration of the expression which is the 
product of the normal density and distribution function in the appropriate 
range of the x values. Specifically, for the computation of F(s), the 


10 


0.50 


Conditional density 
Standard normal density 



00 -2.00 -1 .00 0.00 I'.oo z'.oo 

X 


Figure 3 . Conditional density of X given -«> < Y < t (t = 1 . 00 , 
probability = 0.8413) where (X, Y) is bivariate normal with p = 0.6000, 
and standard normal density . 


value of double integrals such as 


G.(s,t) =f exp{ -x2/2) f 

K Jf-00 4/-0 


(t - px)//l - 


exp( -u^/2)du 


dx 


for all values of s, t and p are required. In terms of these functions, 
it is easily seen that 

F(s) = |2Tr[$(b) - $(a)]| ^ |^Gp(s,b) - Gp(s,a)J 


A closed analytical expression for Gp(s, t) is not available and for 

specific values, numerical methods may be employed. However, in cases where 
symmetry occurs, the numerical computations for a smaller range of values 
are needed. In order to calculate F(s) for all values of s, a, b and p, 
a computer program using quadratures was developed at DFRC and is given in 
the Appendix. 

Rectangle Probabilities. The region (c < x < d, a < Y < b) is a rectangle in 
the (x, y) plane. Thus the joint probability Pr(c <X<d, a<Y<b) 
for real values of a, b, c and d corresponds to a rectangle probability. 

The appended computer program can be used to calculate all such rectangle 
probabilities. The procedure is to identify first that 

Pr[c < X < d, a < Y < b] = Pr[c < X < d|a < Y < b] Pr[a < Y < b] 

= CF(d) - F(c)] Pr[a < Y < b] 

= [F(d) - F(c)] [$(b) - $(a)] 


for all values of c < d and a < b, and then use the computer program with 
the proper inputs. 


COMPUTER PROGRAM INPUTS AND OUTPUTS 

The computer program developed at DFRC computes the conditional density 
and distribution function as outputs for specified values of x given the 
end points of the interval of the conditioning variable Y, and the corre- 
lation coefficient p. Thus the inputs to the program are specific x 
values, end points of the Y interval and the p value. The output has 
two options. Either the density or distribution function, or both may be 
obtained by stating the options in the program. 

The rectangle probabilities are to be obtained by finding the conditional 
probabilities. The computer program with its options is explained in the 
Appendix. 
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EXAMPLES 


The following examples Illustrate the use of the program and tables 
shown in the Appendix to calculate various probabilities. 

The data for the examples are taken from a Closed Circuit Television 
(CCTV) experiment. In this experiment, two pilots, A and B, landed an air- 
craft with the help of an airborne television camera and video monitor. 

Each pilot made ten (10) touchdowns under visual flight regulations, and 
eighteen (18) touchdowns utilizing the closed circuit television monitor. 

The summary of data from the twenty-eight (28) touchdowns is given in Table 1. 
For this illustration the data parameters are vertical acceleration, 

forward velocity, and lateral acceleration Ay. 


TABLE I. SUMMARY OF 28 TOUCHDOWN DATA OF CCTV EXPERIMENT 


Pilot 

Parameter 
(Units ) 

Mean 

y 

S.D. 

a 

Correlation Between 

A 

A^(G) 

1.313 

.2021 

(A^.V^,) = .2481 


v^(mph) 

60.25 

1.3089 

(A^,Ay) = -.0715 


A (G) 
y 

.023 

.1227 

(V^,Ay) = .2807 

B 

A^(G) 

1.294 

.1044 

! 

(A ,V ) = -.2569 
z c 


V^,(MPH) 

62.04 

1.8747 

(A^,Ay) = -.2199 


A^(G) 

-.007 

.0801 

(V^,Ay) = -.1993 


The variables (A , V_, A ) are assumed to follow a multivariate normal dis- 
z c y 

tribution. Thus any two variables follow a bivariate normal distribution 
and any single variable, a univariate normal distribution, as shown in 
figure 1. Further, all the values in these date are considered to be param- 
eter values. 

Example 1 . Computation of a 95% confidence interval of forward velocity (V^) 
given vertical acceleration (Ay) mean is within ± one standard deviation (a). 
It is desired in this example to determine a 95% confidence interval for air- 
craft forward velocity (V^), in miles per hour, at the point of touchdown, 

given the pilot's average vertical acceleration (A^), in G's, within ± one 
standard deviation. The 95% confidence interval end points for V^. given A^ 
mean is within ± a are obtained by solving for t from the equation 




.95 = Pr[-t < (Vj, - Wq)/cTj, < t|-l < (A^ - < 1] 


= Pr[-t < X < t|-l < Y < 1] 


and identifying the interval as (-ta^ + ta^, + y^.). 

The solution of the equation for pilot A data of y^. = 60.25, = 1.3089, 

y^ = 1.313, 02 = 0.2021 and correlation (A^, V^) = -.2481, yields the value 
of t = 1.91666. The 95% confidence interval, therefore, becomes 

(57.7413, 62.7587) 

This shows that if in pilot A data, the aircraft's vertical acceleration 
at touchdown is within +1.3 ±0.2 G's, he has a 95% chance of landing the 
aircraft between 58 and 63 MPH. 

For pilot B data, from table 1, the t value computes to be 1.9136. 

Thus the 95% confidence interval is 

(57.7453, 62.7547) 

indicating if pilot B's vertical acceleration data at touchdown is within 
+1.3 ±0.1 G's, he also has a 95% chance of landing the aircraft between 58 
and 63 MPH. 

Example 2 . Computation of the probability that the forward velocity (V^.) 

and (Ay) are both within ±ct of each variable. The probability of Vj, and Ay 

being within ±a of each respective mean is an example of rectangle proba- 
bility. In this example, the probability that simultaneously, V_ and A.,, 

will be within one standard deviation of each variable's respective mean is 
to be computed. 

This rectangle probability can be obtained by finding 

Pr[-1 < (Vc - yc)/^c < 1 . -1 < (Ay - Vy)/ay < 1] 

= Pr[-1 < X < 1 I -1 < Y < 1] Pr[-1 < Y < 1] 

From univariate tables, Pr[-1 < Y < 1] = $(1) - $(-l) = .6826 and is not 
affected by the correlation coefficients. In order to obtain Pr[-1 < X < 1 
I -1 < Y < 1]^ the values of the correlation coefficients are needed. 

The correlation coefficient (V^., Ay) for pilot A data is equal to -0.2807. 

The computer program output, therefore, for this correlation yields 

Pr[y(, - , Vy - + Oy] = .47554 

Thus, for pilot A there is a 48% chance that simultaneously at touchdown, 
the aircraft's forward velocity will be within 60 ±1.3 MPH and the lateral 
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acceleration is within 0 ±0.1 G's. Conversely, the probability is 0.52 
that both variables will not simultaneously be within one standard deviation 
of their respective means. Similarly, for pilot B with the correlation 
(V^, Ay) equal to .1993, the program yields 

Pr[yc - <^c < Vc < Pc ^c» ^y “ '^y ^ ^y ^ ^y ®y^ " .47078 

which represents a 0.47 probability that the forward velocity will be within 
62 ±1.9 MPH and lateral acceleration is within 0 ±0.08 G's. 

Example 3 . Computation o f the probability of forv/ard velocity (V^) and 

lateral acceleration (A^) being within ±a of each variable, given vertical 

acceleration is equal to its mean (Ay = ^ 7 ) . This rectangle probability can 

be obtained as in Example 2, except in this case the vertical acceleration 
(A^) is set equal to the variable's mean value (v 2 )* probability in other 

words, is a function of a conditional correlation coefficient which is differ- 
ent from the coefficient given in the table. 

For the pilot A data, this conditional coefficient is equal to .3809 and 
the program output yields 

^ '^c ^ ^c °c ’ ^y ■ ''y ^ \ ^ ^y ‘^'y^ = 

which represents a 0.48 probability that the forward velocity will be within 
60.25 ±1.309 MPH, and lateral acceleration within .023 ±0.1227 G's given 
that vertical acceleration is 1.313 G's. 

For pilot B, the conditional correlation coefficient is equal to -.2807 
and the corresponding rectangle probability is 

Pr[yc - < Pc > viy - cjy < Ay < Py + Oy] = .47554 

This represents a 0.48 probability that the forward velocity will be within 
62.0 ±1.9 MPH and lateral acceleration is within 0 ±0.08 G's given the ver- 
tical acceleration is +1.294 G's. 


Dryden Flight Research Center 
National Aeronautics and Space Administration 
Edwards, California, March 17, 1980 
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APPENDIX 


The program to compute the conditional density and distribution 
function for specified values of x given the conditioning on 
variable Y. 


PROGRAM MAIM 


73/74 QPT=1 


FTN 4.2+75360 


5 


10 


15 


Z**** 

z**** 

z**** 

c**** 

C + + + * 


C + + ** 


1 


PROGRAM MAIM(OUTPUT) 

ILLUSTRATIVE USE OF THE ENCLOSED COMPUTER 
PROGRAMS TO COMPUTE VARIOUS 
PROBABILITIES ASSOCIATED KITH THE EXAMPLES 
GIVEN IN THE TEXT OF THIS PAPER... 

BY 3R0.(NL0W^ SOC/ IS I 

PRINT 1. RPR0B(-1.»1.»-1. »!.»-. 2307) 

PRINT 1> RPR0B(-1.>1.»-1.>1.>-.1993) 

PRINT 1> RPR08(-1.»1.>-1.>1.» .3809) 

PRINTl. TINV( . 95.-1. >!.»-. 2401) 

PRINT 1. TINV( .95>-l.>l.»-.2569) 

FORMAT!* *F10.5) 

END 
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TIQM CO 


73/7^ OPT=l 


FTN H. 2+75060 


10/25/79 


+ + + 

ct* + + 
C + + + + 

c*** + 




FJ.^CTnS C0(X) 

COMOITIONAL OEr^SITY FUr^CTIOM 
CD(X\A<y<B) = i./SQRT(2^PI*EXP(X+X))« 

( ^HI( (3-R*X >/ (SORT( 1-R+R) ) - PHK ( A-Rt X ) / ( SQR T ( 1-R+ R ) ) ) 
/( P4K3) - PHKA) ) 

WHERE R = COEFFICIENT OF CORRELATION BETWEEN 
X AND Y 


T 

PHKT) = INTEGRAL F(X) DX 
-INF 

ANO F(X) » l./SQRT(2+PI^EXP(X+X) ) 
BY BROWNLOWf SDC/ISI 


COMNON/PARAN/A, B,R>SQR 

CO = #3989^r228/SQRT(EXP(X + X)) ♦ ( PHK (3-R^X)/SQR)- 
PHK ( A-R*X)/SOR) )/(PHKfl)-PHKA) ) 

RETURN 

ENO 
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FUMCTIGN FI'>IT2 


73/7^ 


o?r = i 


FT*^ 4.2+75060 


10 


15 


20 


25 


30 


35 


40 


45 


50 


55 


Fj:MCrt3>^ FINT2(F»Af3) 

IMTtO^AL OF THE FJHC TION F FROH A TO 9 
+ BY 34J5S I AS-Lif;^FMDRE OLIAOR ITUR E .. 96 POINT FORH 

RtOOlRES 96 E7 ALJ AT IONS QF F(X). 

F must 9E DECLARED EXTERNAL IN 
THE CALLING PR10RAN. 

3Y BROWNLQWf SOC/ISI 

C ^ 

03USLE PRECISIJN POOT ( '13 ) > WEIGHT(43)» AnI S W ER » 0 A, DB 
ODOBLE PRECISION ARY(^3.2) 

EQUIVALENCE ( ARY ( 1 . 1 ). ROOT ( 1 )) > ( ARY ( 1 , 2 ) , RE 1 3HT ( 1 ) ) 

(^ ♦ ^ + ^ 

^ t t 

Sf^T UP ROOFS ANO W EIGHTS, .. 

•DATA ( ( ARY( I . J ) * J =1» 2 ) » 1 = 1 » 16) / 

. 3.31327 o7AA3 A9602 36Vj79aO> 0.03255 0614A 923b3 1662A2D0. 

. O.OA331 29351 360A9 73111200,. 3.33251 51187 13868 335937D0, 

. 0.05129 7A9bA 5'-tA2p :>y899A00> 0.032A4 71637 14064 269364D0. 

. 0.11369 53501 10565 9209110 0. 0.03234 38225 66575 92342900. 

. 0.14597 37146 54396 94198930. 0.03220 62047 94030 250o69D0.. 

. 0.17309 63823 57513 60275900, 0.03203 44562 31992 66321300, 

. 0.21003 13104 60567 20360300, 0.03132 37533 94411 00653500, 

. 0.24174 31361 63343 31232300, 0.03158 93307 70727 160558DO. 

. 0.27319 33125 91349 14143700, 0,03131 64255 96361 35531300, 

. 0.30436 49443 5‘.496 35302403. 0.03101 03325 85313 63742300, 

. 0.33520 35223 92625 42261600, 0.03067 13761 23669 14901400, 

. 0.355o9 68614 72313 63533103, 0.03029 99154 20327 59379400, 

. 0.39579 76493 28903 603235J0, 0.02939 63441 36326 38593400. 

. 0.42347 39334 37300 34536500. 0.32946 10399 53167 90597000. 

. 0.43t70 94221 67743 00863603. 0.32399 46141 50555 23654300, 

. 0.43345 79739 20596 359768D0, 0.02849 74110 65085 38564600, 

. 0.51169 41771 54667 67353600, 0.02797 00076 16348 334440DO, 

. 0.53933 31033 24357 43622700, 0.02741 29627 26029 24232300/ 

C ^ 

OATA ((ARY(I,J),J=1, 2). 1=19,37)/ 

. 0.56661 04135 61397 16340400, 0.02632 68667 25591 76219300, 

. 0.59303 23647 77572 03363400, 0.02621 23407 35672 41391300, 

. 0.61392 58401 25468 57033600, 0,02557 00360 05349 36149900, 

. 0. 04416 34037 34967 10679300. 0.02490 06332 22483 61023300, 

. 0.66371 33100 43916 15395300, 0.02420 43417 923o4 69123200, 

. 0.69256 45360 42171 5ol34400, 0.02343 33990 35926 219342D0, 

. 0.71567 03123 43967 62622500, 0.02273 70696 58329 37400100, 

. 0.73303 0643 7 44‘.00 1 32351DO, 0.02196 6o444 33744 34919500, 

. 0.759o0 23411 76647 49870300, 0.02117 29398 92191 29398300, 

. 0.78036 90433 67433 21760400, 0.02035 67971 54333 32459500. 

. O.SOOiO 37441 39140 31722900, 0.01951 90311 40145 02241000, 

. 0.31940 03107 37931 67553900, 0.01366 06796 27411 46733500, 

. 0,337o2 35112 23127 12149400, 0.01778 25023 16045 26083390, 

. 0.35493 70334 34601 45546300, 0.01683 54798 64245 17245000, 

. 0.37133 35059 09296 502B74D0, 0.01597 05629 O25o2 29138100, 

. 0.33639 45174 02420 41605700, 0.01503 87210 26994 93300600, 

. 0.90146 06353 15852 34131700, 0.01409 09417 72314 36 3 91600, 

. 0.91507 14231 20393 07420600, 0,01312 32295 6o961 57263700. 

. 0.92771 24567 22303 69096500, 0.01215 16046 7103B 31963500/ 

CT' + t'Y 


19 



FUNCriG‘>i FInIT 2 73/74 GPT = l PIN 4.2 + 7t?060 




OArA( (ARy(i^j)>j 

= 1.2 ) • 

1=33,43)/ 







0.93937 

03397 

52 755 

21693200, 

0.01116 

21020 

99338 

493591D0. 

60 


0.95003 

27177 

34937 

63575600. 

0.01016 

07705 

35003 

41573800, 



0.9396d 

32914 

96742 

53930090, 

0.00914 

36712 

30733 

38663300, 



0.96332 

632B4 

63264 

21217400. 

0.00312 

63769 

25698 

75921700, 



0 . 9759 3 

:^L745 

35136 

46645300, 

0.00709 

64707 

91153 

86526900, 



0.93251 

72635 

63014 

67744700. 

0.00605 

35455 

04235 

96168300, 

65 


0 . 9 3 d 0 5 

^1263 

2962 3 

79948100, 

0.00501 

42027 

42927 

51769300, 



0.99254 

39003 

23762 

62^57200, 

0.00396 

45543 

33444 

68667400, 



0.9 ?593 

13429 

37209 

29065000, 

0.00291 

07318 

17934 

94640300. 



0.9^336 

43753 

o3161 

67772403, 

0.00185 

39607 

38946 

92173200, 



0.99963 

95 03 b 

33230 

76662800, 

0.00079 

67920 

65552 

01242900/ 

70 





















c ^4^ t > ^ 


INTEGRATIJN DONE BY 

TRANSLATING F 

TO THE 





INTERVAL 

-1 TO 1 






C « ^ « 










75 = 0.90 

93 * 3 
OA = A 

03 i I=l>43 

83 r * { ( 93-OA) ♦ROOT ( I ) + ( 0 B+ D A ) ) / 2 . DO 

A^iS^ER = A>15>JER + ^£I3HT(I) + F(T) 
r = ( (03-DA) + (-R3GT( I) ) + { 9 8 ^- D A ) ) / 2 . DO 
1 A454ER = A'^SWER + y^EIGHT(I) + F{T) 

C =4' 

85 FIMT2 = {9^-aA)+AMSWER/2.D0 

C 

RETJRM 

£nIO 
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FUNCTIGM RECT 


73/7^t 


GPT = 1 


FTN 2+73060 


+ + 

+ + * 


C + + ** 


C + + ** 

c + +++ 


FU'^CTION fttCT(A,a.rO 


RFCTA'IGLE P <"33 A 3 1 L I T Y . . . 

V^LU'IE 'JSOER the NTRHAL BIVARIATE DENSITY* 
-INF<X<Af -INF<Y<3. 

BY 3R3WNL0W* SOC/ISI 


C0HM0N/6PAR>1/ AA,B3*RR>SaR 
EXTERNAL 6 

A A * A 
36 = 3 
R R = R 

SOR » SORT(i,-R+R) 

RECT = FIMT2(G>-13.f A) 

RE TURN 
EHO 



FUNCTION G 


73/74 GPT=1 


FTN 4.24^75060 


FUNCTION G(X) 



C + + + * 

conditional distribution FUNCTION... 


C + + + + 


5 


COiMON/GPARN/A,B,R,SOR 

10 

C** + 4e 

T = (3-R*X)/S0R 

G * EXP(-X*X/2. )*?HI (T)*2. 506628275 

RETURN 

END 


22 



FUMCriGM TIN7 73/7^ jPT=l 


FTM ^,2+75060 


1 


10 


15 


20 


25 


30 


35 


FJ.'JCTinN TINVt^f A,Q>R) 

GIVEN A<Y<6 FIND T SO THAT -T<X<T AND 
:♦♦♦♦ ?( -T<X<T \ A<r<3) = P 

WITH COEFFICIENT OF CORRELATION BETWEEN X AND 
Y cDUAL TO R« 

C * =45 ^ + 

+ + P( -T<X<T \ A<Y<B) = P( -T<X<T» A< Y< B ) / P ( A < Y < B ) 

c ♦ ^ + 

T IS FOUND 3Y INTERVAL HALVING. 

IROWNLOWf SDC/ISI 

DENON = PHI(3)-PHI(A) 

TiAX » 10. 

TilN = 0. 



■J ] 1 1 = 1. 5D 





r » (TNAX + TNIN)/2. 




PCJiPT = RPR03(A. 

3r-T,T 

,R) /OE 

NGN 

C > 4c 

IFIPCONPT .&r. P) 

TNAX 

= T 



IFCPCO'IPI .LT. ?) 

TM.IN 

= T 



IF( A3S( PCOMPT-P) 

.LE. 

l.E-5) 

GG TO 2 


1 CJNTINUE 

PRINT 100, P.A,3,R 

100 F.)RMAT(* COULDN’T FIND T IN 50 ITERATIONS: P»*>F7.4,+ A»*, 

, + B=*,F7,', Ra*tF7«A) 

TINV = { TNIN^TMAX ) /2. 

RETURN 

C + + 

2 TINV = T 
RETURN 
END 


F7.4t 
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FJ‘4CTIJM 


73/7^ 




PTN 4.2+75060 


5 


10 


15 


^eCTAlbLz j^=!13A3rLlTY FO-^ BIVA-^IATE 
:♦♦♦♦ NJR>!AL OISr^IOJTION... 

C + + + ^ 

C + C<X<0 

A<Y<3 

c + *t+ AND T^F COEFFICr-NT TF C L A T I GN BETWEEN 

X AO Y IS 

3r i^DWNLJW. S0C/I5I 

!^Pr^G3 = (RiCnO.Bf-?) - RECTC,3.R)-^ECT(D^A,R) + RECT(C>A#R)) 
. ^0.15^l5t045 

c ♦ t ♦♦ 

RETURN 

END 
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FUNCTI3M PHI 


73/74 


OPT* 1 


FTN 4.2+75D60 


RETURN 

tSD 


FJNCTIO'i PHI 


73/7^ 


J?T = 1 


fJH ^.2+75060 


5 


10 


15 


20 


25 


30 


35 


^0 


^5 


50 


55 


FjNCTiJM pnnx) 

NGRMALOfl) DISTi^IBUTiaN FUMCTION 
PHKX) = rUEGRAL GF NORMAL DE!>ISITY 
FROM -INFINITY TG X. 

C + + ** GY -jRJvJNLJWf SOC/ISI 

r ♦♦♦♦ 

C^^ + + 

LOGICAL FLAG 

IF (X .GT. -iO, ) GG TG 1 

PHI » 0. 

RETURN 

+ + + 


1 


c«* + + 

2 

C + t* + 
C + =4^ + + 

C ♦ + ♦* 

c + + t + 


z^*** 

3 




C + =4= + + 

c + + ^ + 

C X 


C ♦* 
C + + ♦ + 


C + 

C ’J' * ♦ -Jf 

c 

C + + + + 

C’4^ + * + 
5 

C + + + + 


IF(X .LT. 10.) G3 TG 2 
PHI = 1. 

RETURN 

FLAG = .T. 

JETcRMlNE IF X>0» SERIES EXPANSION IS FOR 
POSITIVE VALUES GF X.. 

IF(X .GT. 0.) GG TG 3 
FLAG = . F . 

INITIALIZE VALUES FOR PARTIAL SUM OF THE SERIES... 

Z = A3S( X) 

0 = 1 . 

SUM = 0. 

TJP » Z 
3JT = i. 

CGNTI NUt 

SAVE = SUM 

SUM » SUM 1 - rOP/6GT 

CONTINUE TG SUM UNTIL MACHINE UNDERFLOWS... 

IF(SAVE .E3. SUM) GO TO 5 

UPDATE EXPRESSIONS FOR THE SUM... 


TOP = TG?+Z«Z 

0 = 0 + 2 . 

3JT = 3GT+G 
G3 TG ^ 


DEPENOING JPQN WHETHER ORIGINAL X>0 OR X<0> 
GET APPROPRIATE INTEGRAL VALUE... 

•^HI = SJM/5QRT{6.2331553G3 + EXP( X + X) ) +.5 
IF(FlAG) return 


=>HI = I.-PHI 
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FUNCTION FRA 


73 / 7 '^ 


GPT = 1 


FTN 4.2+75060 


( 


5 


10 


15 


20 


25 


30 


35 


40 


45 


FUNCTIOiM FRAC(PV, Yi,Y2) 

C + + + + 

C+++* GIVEN A 3IVARIATE NORMAL OISTRISUTION 

^ITH COEFFICIENT OF CORRELATION RHO 
C + * + + AND Yi < Y < Y2» F R AC < V^ Y 1. Y2 ) RETURNS 

C+*** THAT VALUE T# SUCH THAT: 

C ♦ + + ’4' 

C++** PR03( -T < X < T, Yl < y < Y2 ) = PV 

C + + + + 

C++++ BINARY SEARCH^ LIMITED TO A MAXIMUM OF 20 ITERATIONS 

C + + + + 

C + + + * 

C++++ BY 3R3WNL0W> SDC/ISI# 11/79 

C + + + + 

C + + + + 

KOUNT = 0 
TMIN * 0. 

TMAX = 10. 

C + + ** 

1 T * ( IMIN + TMAX) /2. 

C + + + + 

VAL * RtCT(-T,T#Ylf Y2) 

PRINT lOOf VAL*T 
100 FORMAT(+ +2F10.6//) 

C + + + + 

C + + ** 

C++++ IF WERE WITHIN l.E-5 OF THE VALUE# WE 

C+++* HAVE FOUND THE SOLUTION... 

C + + + + 

IF( ABS(VAL-PV) .LT. l.E-5) GO TO 2 

C + + + + 

IF(VAL .LT. PV) TMIN=T 
IF(VAL .GT. PV) TMAX « T 

C + + + + 

C++** CHECK FOR MAXIiUM NUMBER OF ITERATIONS... 

C + + + + 

IF(KOUNT .GE.20) RETURN 
K3UNI * KOJnT + 1 

C + + + + 

c + + + + 

GJ TO I 

C + + + + 

2 FRAC * T 
C + + + * 

C + + + + 

RETURnI 

E>IO ‘ ' 
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FUNCTION C3NDEN 


7i/7^ 


GPT = i 


FTN 4.2+75060 


C + t + t 

+ + ^ 
C + + + + 

c ** + ’»' 

c + ++« 

C + + + + 
C + + + + 


c^ * + ♦ 


function C‘JNDEN()() 

CONDITIONAL DENSITY FUNCTION OF X, GIVEN 
A<Y<S FRON BIVARIATE NORMAL DISTRIBUTION 
F(X,Y) WITH COEFFICIENT OF CORRELATION RHO. 

BY BROWNLOW^ SOC/ISIf 11/79 


CDMMON/ARG/RHG, A#8 

SET UP THE PARAMETERS^ PHI IS THE UNIVARIATE 
NORMAL DISTRIBUTION FUNCTION. 

R - SORT(1.-RHO*RHO) 

D * PHI(B)-PHI ( A) 

T = PHK (3-RH0tX)/R) - PHK (A-RHO*X)/R ) 

CONDEN = EX?(-X^X/2. )tT/(D*2. 506628275) 

RETURN 

END 


FUNCTIOM 3 


73 /7i 


J ^ T = I 


FTN ^.a+75060 


C + ^ 
c + + + t 

C¥¥*^ 

C+ 

c « + > « 
c« + + * 

C>4'*4' + 

c« +=!« + 
ct =«'♦♦ 
C * 

C + + + 4' 
c + ♦ t ♦ 


r J 4 C T i InI G ( 5 I T > 

AWA'iUTe .'I3=,'IAL DIST»I3UTinS FUNCTION. 
O(S.l) = OQJHLt INTEGRAL OF NORMAL 
>)Ii/AKlArt OENilTf FUNCTION, -INF TO 3, 
-INF TO T. 

NJTUE THAT THE NUMERICAL COMPUTATIONS 
U3c THE FACT THAT THE CONTRIBUTION TO THE 

integral nalur from -inf to -15. 

13 INSIGNIFICANT. 

AT HRONNLOR, SOC/ISI, 11/79 


COMMON/RAjS/IT 
EXTERNAL F;N 
TT = T 

G = FINT2(FCN.-15.,S)/3.?33135308 


return 

ENJ 


FUNCTIQM FCU 


7i/7^ 


■]PT = 1 


fJH 4.2+75060 


FJNCTIJN F:N(X) 




O^NSirv FUMCTIG:^ FGR DJJ9LE INTEGRAL, 


C t 

DHI(X)+Z(X)f WHERE PHI AND Z ARE THE 

5 

C + + + + 

NjK'iAL 01STR13'JTIQN AND DENSITY FUNCTIONS 


C + + ** 

RESPECTI VELY. 


C + + ** 

C^t*** 
C + + ** 

^Y iPl^NLD^, SUC/ISI, 11/79 

10 

c + ++ + 


15 

C + + ** 

C + + ** 

CJifllN/ARG/PHOf 
C 1NNON/PA55/TT 

Z(ARG) = [ XP(-aP3+AP3/2. ) 

20 

C + + ** 
C+ + ** 

C + + ** 

0 = ( TT-RH J+X ) /SOP F( I.-RH3+9HG) 
FCM = PHK J)+Z(x) + 2.:iOD62S275 


C + + +♦ 

c ♦ + + ♦ 

RETURN 

c'iO 


30 


FJNCTn^) ^-.CT 


73/7'^ 


G^T = i 


FTN 4.2+75060 


0 


5 


10 


C + + t * 

C + + + + 
c + + * ♦ 

C + + *4t 
Q + + « « 

C ’i' + * 

C + **1C 


cjnCTIJN «!ECr(Xl»X2,n»Y2) 

'^eCTi'liLc P?J3AnulTIES FOP BIVARIATE 
DlSTKiBJTIO'^. <1<X<X2* Yl<r<Y2> AND 
CTEt-FICIE'IT JP CORRELATION IS RHO. 

3Y BRl^fNLOR* 3DC/ISI» 11/7R 


RECT* G(X2.Y2»-5(<i»Y2t-G(X’.Yi)+G(Xl*Yl) 

RETURN 

END 


NORMAL 

THE 


31 



rUNCTluN CONJIST 




JPT = 1 


FTN 4.2^*75060 


0 


5 


10 


15 


FJ.^CTIOM C1^>10IST ( <) 

C + + + + 
c« + 

C J.l.i} J/AKt/PHJf A.r^ 

PHI 13 THE U.SIVA^IATE NOkHAL 0 I S TR I U T I QN 
+ FUNCTiQNl. 

C + + 

C0H015 r= (G( X, H)-, (X» A) > / { (PHI (3)-PHI ( A ) ) + 6.2d318 53 08) 

PETURH 

FHO 


CaHOITIOHAL DISTRIBUTION FUNCTION OF X GIVEN 
A<Y<3 FROi BIVARIATE NORMAL DISTRIBUTION 
FCXfY) WITH COEFFICIENT OF CORRELATION RHO. 

BY 3R:lwNL0W> SOC/ISI, 11/79 
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FUNCTIG'^ FIST? 


7 3 / 7 ^ 


.1 ^ T = I 


FTN ^.2+75060 


5 


10 


15 


20 


25 


30 


35 


^0 


45 


50 


55 


FJ'JCTIJ'I FlNr.2(F,A*3) 

C ^ 

OF THF FUNCTION F FRQI A TO B 
-M r,.VJj >I AN-L£">eNO?E QU40R1T'JRE> ?6 POINT FIRM 
C + 1'** REjUIRCJ TE evaluations of F(x), 

c+ ♦ ♦ * 

C«*M* F “iUST BE OF.CLARC-D EXTERNAL IN 

C*»*+ TH-r CALLUS 0R.K-.RAM, 

CM*** BY IF Ow 'ML OR* SOC/ISI 

CM*** 

OIJMLc PRtCIjlOi RJOT(Ad)» VEI3HT(A3)» ANS 'W £ R » 0 A» 0 B» AR Y ( 4 3 » 2 ) 
EJUIVALE'MCE (ARY(l»i)>ROOT(l) )» ( AR Y ( I ,2 ) . RE 13 MT ( 1 ) ) 

C**M* 

c->*** 

C**M* SET JP R-J0T5 AND HEIGHTS... 

UATA ( ( AR Y ( 1 . J ) • J = 1 . 2 ) • I =l * 16 ) / 

. 0.01O27 37AA-J AV602 Y6-757700. 0,03255 061AA J2363 16624200f 

, O.OA'Jji 1V351 360A? 73U120U» U. 03251 61137 13366 335V8700, 

, 0.0312V 7+'V56 6A'i25 dSsOVAOOj 0.032AA 71637 1A066 26V36A00, 

. 0.1136y 'j350, lu665 620911i)J» 0.03236 36225 68575 92362VD0> 

. j.x65'?7 3 7166 56696 961V89D0. 0.03220 62067 96030 250669D0» 

. 0.17607 66823 67613 60275000, 0.03203 66562 31992 663213D0, 

. 0.21003 13106 60567 20360300. 0.03162 37533 96611 006535D0, 

. 0.261/6 31561 63660 01232600, 0.03158 93307 70727 163553D0, 

. 0.27319 33125 91069 16163700, 0.03131 66255 96861 355313D0, 

. 0.306 39 .,9=.93 69996 3p302600, 0.03101 03325 36313 837623D0, 

. 0.33p20 6522m 9?o26 92261o00, 0.03067 13791 23669 16901600, 

. O.J65j9 68616 72313 63503100, 0.03029 99156 20327 593796D0, 

. U.3957.X 76698 2m906 60323500, 0.02989 63661 36328 33598600, 

. 0.62567 89686 07600 p6p36600, 0.02966 10399 58167 90597000, 

. 0.65670 96221 67763 008b3600, 0.02899 66161 50555 23656300, 

. •0.683.,: 79739 20 396 35 976300, 0.02869 76110 65085 38566600, 

. 0.3116 9 61771 t>-,667 6 7 35860 0, 0.0279 7 00076 16868 3366 80 00, 

. O.P393-3 810-3 3 28367 63d22 700, 0.02761 29627 26029 26282300/ 

C**M* 

OATA {(A8Y(1,J),J=1,2),I=19,3Y)/ 

. 0.5o66I 06186 61397 16860600, 0.02662 63667 25591 76219300, 

. 0.69309 23667 77p72 0.j363600, 0.02621 23607 35672 61391300, 

. 0.61-i9'’ 5 )601 2 6 665 57033600, 0.02557 00360 05369 3616 9 9D0, 

. 0.6661m 3-.)37 399m7 1 0679800, 0.02690 06332 22683 61028800, 

. 0.6p87i 83100 6391m 15395303, 0.02620 63617 92366 69128200, 

. '1.5928-p 9,366 -.2171 56 1 39900, 0.02369 33990 35926 21986200, 

. 0.71567 63123 69767 32622500, 0.02273 70696 58329 37600100, 

. ■0. 7 38 J3 06637 69 , Oj 1 328 51 03, 0.021 96 66666 33766 36919500, 

. 0.75760 2 861 1 7669 / 6987030.0, 0.02117 29390 92191 29S98300, 

. O.IiJii) 90938 67,33 21 760900, 3.02035 67971 56333 326595D0, 

. 0. 80030 87661 37160 81722900, 0.01951 90311 60165 02261000, 

. '0.8 1 960 031.J 7 3 7 781 67553900, 0.01366 06796 27611 66738500, 

. 0.837j2 35112 28127 12167900, 0.01776 25023 16065 26083300, 

. 9.35698 90336 36601 95566300, 0.01688 56798 66265 17265000, 

. 0.3713M 85057 09?9o 50237600, 0.01597 05629 02562 29138100, 

. 0.0868 9 9-5 1 79 0?920 6160-5700. 0.01503 8 7210 26996 938006D0, 

. 0.90166 05353 1p352 36131300, 0.01609 09617 72316 36091600, 

. 0.9iP07 16231 20, -7 3 076206DO, 0.01312 82295 66961 57263700, 

. 3.92771 26,67 22303 59096p00, 0.01215 16066 71088 31963500/ 

C**M* 

0AIA((ARY(I,J),J»1,2),1«33,68)/ 
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FU'^CTION FI'IT? 


73/7i. 




FTN A. 2+75060 


60 


63 


70 


75 


80 


85 


. 0.V373? 033';i7 j2753 21‘5932D)» 0.01116 21020 99838 49859100. 

. 0.-/jo0 3 27177 3443 7 o3375600. 0,01016 07705 35008 41575900. 

. 0 . 954,33 32 91 4 48 742 53930000. 0.00914 8o7l2 30783 38663300, 

. ^.9o 3J2 63234 332o4 21217400. 0.00812 68769 25698 75921700, 

. 0.9 73,3 9 i 74 5 35 133 46o4530'J, 0.00709 64707 91153 36526900, 

. 'J.98’31 72 335 53014 577447J0. 0.00605 85455 04235 96168300, 

. 0.-,3305 41 263 2932 3 79943100, 0.00501 42027 42927 51769300. 

. 0.9^23', 3 9 00 3 2 3 ^62 62437200, 0.00396 45543 38444 68667400. 

. 0.VJ593 0 7209 29055000, 0.00291 07318 17934 94640800, 

. 0.99833 43753 oJl'U 67772400, 0,00135 39607 33946 92173200, 

. 0.99938 9-J03 3 83230 76682800, 0.00079 67920 65552 01242900/ 

C* + ♦♦ 

c + * + ♦ 

Z*-*** 1 NlfFG-lAriO-l DOME 8 Y TRAiNlSL ATIN3 F TO THE 

iNIF^/tL -1 TO 1 

C + *1‘4‘ 

Ali^r.2 * C.JO 
03 = 3 
3\ = A 

C9-l‘ + + 

03 1 1=1,43 

T = t ( 98-OA ) ♦P )3T t I ) +(D3 + 0A) )/2.00 
A-l3fciK = . vEIG+TtI) * P(T) 

r = ( ( ■)3-0 A) ♦ (-9 j j r ( I ) ) + ( D8+0A ) ) /2. DO 
1 A'1S.(-k' = A<i3vc 9 + VtlvjHTd) 9 F(T) 

Fi 1T2 = ( )8-0rt) + A>l6Wfc9/2.D0 

C ♦ ♦ + 

8 E T JS\i 
t 30 
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FUNCTIGM 0?T = 1 FTN 4.2+75060 


FUMCTIOn PHI(X) 

C + + ** 

C + + ** 

C + + ** N0P'1AL(0»1) DISTRIBUTION FUNCTION 

5 C++++ PHI(X) . INTEGRAL OF NORNAL DENSITY 

FRO'i -INFINITY TO t. 

LOGICAL FLAG 

10 IF(X .GT, -10. ) GD TO 1 

PHI « 0. 

RETURN 

C + + + + 

1 1F(X *LT. 10.) GO TO 2 

15 PHI * 1. 

RETURN 

2 FLAG * .T. 

C* + 

20 C**=«'* OETERMINE IF X>Of SERIES EXf>ANSIQN IS FOR 

PQSITIYE VALUES OF X.. 

C* + + * 

IF(X .GT. 0. ) GO TO 3 
FLAG = .F. 

25 

INITIALIZE VALUES FOR PARTIAL SUM OF THE SERIES... 

C + ^ + % 

3 Z « A3S(X) 

0 = 1 . 

30 SUM = 0. 

TOP « Z 
30T = 1. 

^ CONTINUE 

35 SAVE = SUM 

SUM » SUM TOP/OGT 

CONTINUE TO SUM UNTIL MACHINE UNDERFLOWS... 

^0 IF(SA;E .hO. SUM) go to 5 

C+++* UPDATE EXPRESSIONS FOR THE SUM... 

C + « + * 

TOP = TJP^^'Z^Z 

^5 0 » D + 2. 

30T = JOT^O 
GO TO ^ 

50 DEPENDING UPON WHETHER ORIGINAL X>0 OR X<0> 

GET APPROPRIATE INTEGRAL VALUE... 

+ + * 

5 PHI = SUM/SQRT( 6.2631d5308*EXP(X+X) ) +.5 

IF(FLAG) RETURN 
55 

°NI = l.-PHl 


35 


FUNCTIO>l PHI 7i/7<» OPT = i 


60 




RETURN 


FTH 


.2+75060 C 
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